##Load packages
library(ordinal)
library(haven)

##Load data
datos <- read_sav("TCR_efficacy.sav")
datos2<-as.data.frame(datos)

datos2$imple_level<-factor(datos2$imple_level, ordered = TRUE)
datos2$country<-factor(datos2$country)
datos2$TT_Defeat<-factor(datos2$TT_Defeat)
datos2$Executive_power<-factor(datos2$Executive_power)
datos2$inst<-factor(datos2$inst)
datos2$leg<-factor(datos2$leg)
datos2$const<-factor(datos2$const)
datos2$crim<-factor(datos2$crim)
datos2$ind_rep<-factor(datos2$ind_rep)
datos2$col_rep<-factor(datos2$col_rep)
datos2$sym_rep<-factor(datos2$sym_rep)
datos2$mat_rep<-factor(datos2$mat_rep)
datos2$non_rep<-factor(datos2$non_rep)
datos2$edu<-factor(datos2$edu)
datos2$other<-factor(datos2$other)


##Multilevel multinomial logistic regression

resmultinew<-clmm(imple_level ~(1|country) + 
                 number_recs +
                 Executive_power + 
                 TT_Defeat + 
                  Polity_score + 
                   CS_Consultation +
                     inst + leg + const + crim +
                    ind_rep + col_rep + sym_rep + mat_rep +
                      non_rep +
                   edu + other, data=datos2)
summary(resmultinew)
oddsrat<-exp(confint(resmultinew))
oddsrat
tablafinal<-cbind(oddsrat,exp(resmultinew$coefficients), summary(resmultinew)$coef)
tablafinal


##Plot confidence intervals
ic_ord<-tablafinal[-c(1:2),c(3,1,2)]
minimo<-min(ic_ord)
maximo<-max(ic_ord)
ic_ord<-data.frame(cbind(rownames(ic_ord),ic_ord))
colnames(ic_ord)<-c("ID","exp_coef","L","U")

ic_ord$ID<-c("Numer of recommendations",
             "Executive Power-Support","Tanstition Type-Defeat",
             "Polity Score","Civil Society Consultation",
             "Institutional","Legislative","Constitutional","Criminal",
              "Individual reparation","Collective reparation","Symbolic reparation",
              "Material reparation","Non repetition","Educational","Other")
          

ic_ord$exp_coef<-as.numeric(ic_ord$exp_coef)
ic_ord$L<-as.numeric(ic_ord$L)
ic_ord$U<-as.numeric(ic_ord$U)
par(las=1)
par(cex=1)
par(mar=c(6, 12, 2, 2))
plot(ic_ord$exp_coef,16:1, main="Confidence Interval", 
     xlab="Odds Ratio",ylab="",pch=17, 
     xlim=c(minimo,maximo), yaxt="n", family="serif")
segments(ic_ord$L,16:1, ic_ord$U,16:1, lwd=1.5)
abline(v=1,lty=3, col="grey")
axis(2,at=1:16, labels=ic_ord$ID[16:1], cex=0.05, las=2,family="serif")
